#include "header.h"

//------------------------------------------------------------//

void dump(FILE* fh, struct grid* in, double* time){
  //writes grid to file

  int j;
  double dd2, dr, d0, df, db, dp2;		
  double phase0, phaseb, phasef;		

  for (j=2; j<(in->Np)+2; j++){ 
    dr = in->dr;	
    d0 = sqrt(pow(in->ugrid[j], 2) + pow(in->vgrid[j], 2));
    df = sqrt(pow(in->ugrid[j+1], 2) + pow(in->vgrid[j+1], 2));
    db = sqrt(pow(in->ugrid[j-1], 2) + pow(in->vgrid[j-1], 2));
    dd2 = (df - 2.*d0 + db)/(dr*dr);
    phase0 = atan( (in->vgrid[j])/(in->ugrid[j]) );
    phasef = atan( (in->vgrid[j+1])/(in->ugrid[j+1]) );
    phaseb = atan( (in->vgrid[j-1])/(in->ugrid[j-1]) );		
    dp2    = (phasef - 2.*phase0 + phaseb)/(dr*dr);	

  		
    fprintf(fh, "\t%.15e\t\t%.15e\t\t%.15e\t\t%.15e\n", in->rgrid[j],d0 ,in->ugrid[j],in->vgrid[j] ) ;
  }  
  fprintf(fh, "\n");
}

//------------------------------------------------------------//

void read_p(FILE* fh, struct INPdata* input){
  //reads all input parameters 

  FILE*   fh2;
  char    line[80], value[80];
  char    runname[80];
  int      n, i = 0;
  
  strcpy(runname, input->runname );
  strcat(runname, ".log" );
  fh2 = fopen( runname , "w" );
  fprintf(fh2 ,"This file contains parameters that were read at execution\n\n");

  while(fgets(line, 80, fh) != NULL){
    i++;
    sscanf(line, "%s", value);
    if (i==2){
      (input->A) = atof(value);
      fprintf(fh2, "A\t =\t %f\n", input->A);
    }
    if (i==3){
      (input->Ngrids) = atoi(value);
      fprintf(fh2, "Ngrids\t =\t %d\n", input->Ngrids);
    }
    if (i==4){
      (input->dr) = atof(value);
      fprintf(fh2, "dr\t =\t %.5f\n", input->dr);
    }
    if (i==5){
      (input->lambda) = atof(value);
      fprintf(fh2, "lambda\t =\t %.5f\n", input->lambda);
    }
    if (i==6){
      (input->tinc) = atof(value);
    }
    if (i==7){
      (input->tmax) = atof(value);
      fprintf(fh2, "Termination time is\t %.5f\n\n", input->tmax);
    }
    if (i==8){
      (input->Gaussian) = atoi(value);
      fprintf(fh2, "Type of IC: \tGaussian\t =\t %d\n", input->Gaussian);
	}	
    if (i==9){
      (input->sigma) = atof(value);
      fprintf(fh2, "Sigma\t =\t %.5e\n", input->sigma);	
	}
    if (i==10){
      (input->Noutfile) = atoi(value);
      fprintf(fh2, "\nWrite selfsimilar data every %d steps\n", input->Noutfile);	
      fprintf(fh2, "Write profiles every %.2f sec\n\n", input->tinc);
	}
    if (i>11){
      n=i-12;
      (input->Npoints[n]) = atoi(value);
      fprintf(fh2, "Grid %d\t\t%d gridpoints \n", n, input->Npoints[n]);
    }
  }
  fclose(fh2);
}

//------------------------------------------------------------//

void write_out(FILE* fh, struct grid* in ){
  //writes to binary file for restarts
  
  fwrite(in->rgrid, sizeof(double), in->Np+4, fh);
  fwrite(in->ugrid, sizeof(double), in->Np+4, fh);
  fwrite(in->vgrid, sizeof(double), in->Np+4, fh);
  //fwrite(in->dgrid, sizeof(double), in->Np+4, fh);
}

//------------------------------------------------------------//

void read_in(FILE* fh, struct grid* in){
  //reads from binary file for restarts

  fread(in->rgrid, sizeof(double), in->Np+4, fh);
  fread(in->ugrid, sizeof(double), in->Np+4, fh);
  fread(in->vgrid, sizeof(double), in->Np+4, fh);
}
//------------------------------------------------------------//

void writeall(FILE* fh, struct grid *m, struct SSdata *S, int *Ngrids){
  struct grid* temp = m;

  fwrite(Ngrids, sizeof(int), 1, fh);

  fwrite(&(S->Phtime), sizeof(double), 1, fh);
  fwrite(&(S->Sstime), sizeof(double), 1, fh);
  fwrite(&(S->L), sizeof(double), 1, fh);
  fwrite(&(S->A), sizeof(double), 1, fh);

  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}

  while (temp != NULL){
    write_out(fh, temp);
    temp = (temp->rpointer);
  }
}
//--------------------------------------------------------------//

void sgaussIC(struct grid* in, int* Ngrids, double* A, double* sigma ){
  /*int j,k,N;
  struct grid* temp = in;
  double i0 = 0;
  double dr;	
  double C0 = 8.95877217377748;
  double fa,fb,fab,f3ab,fa3b;


  for (j = 0; j < *Ngrids; j++) {

    temp->mgrid[2] = i0; 
    N = temp->Np;
    dr = temp->dr;	

    for (k = 2; k < N+2; k++ ) {

      fa =    (temp->rgrid[k])                   *   exp(-pow( (temp->rgrid[k])           / (*sigma), 6 ));

      fb =    (temp->rgrid[k+1])                *   exp(-pow( (temp->rgrid[k+1])        / (*sigma), 6 ));

      fab =   (temp->rgrid[k] + 0.5*dr)        *   exp(-pow( (temp->rgrid[k] + 0.5*dr ) / (*sigma), 6 ));

      f3ab =  (temp->rgrid[k]  + 0.25* dr)     *  exp(-pow( (temp->rgrid[k] + 0.25*dr) / (*sigma), 6 ));

      fa3b =  (temp->rgrid[k]  +  0.75*dr)     *  exp(-pow( (temp->rgrid[k] + 0.75*dr) / (*sigma), 6 ));

      temp->mgrid[k+1] = temp->mgrid[k] + (*A/4)*(1./90)*C0*(dr)* (7*(fa+fb) + 32*(f3ab+fa3b) + 12*fab);

    }
    i0 = temp->mgrid[N+2];
 	
    temp = temp->rpointer;

    if (temp == NULL) break;
  }
  */
}
